Setting the data
# Get palestinian cities, convert it to sf dataframe, assign projection to the coordinates and save it
library(maps) # Get palestinian cities
## Warning: package 'maps' was built under R version 4.4.2
##
## Attaching package: 'maps'
## The following object is masked from 'package:purrr':
##
## map
library(sf) # convert it to sf dataframe, assign projection to the coordinates and save it
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
data(world.cities)
palestine_cities <- world.cities %>% filter(country.etc == "Palestine")
class(palestine_cities)
## [1] "data.frame"
palestine_cities <- as_tibble(palestine_cities)
write_excel_csv(palestine_cities, file = "D:/R_Gaza/pal_map/data/palestinians_cities.csv")
names(palestine_cities) #get the names for the coordinate variables
## [1] "name" "country.etc" "pop" "lat" "long"
## [6] "capital"
palestine_cities_sf <- st_as_sf(palestine_cities, coords = c("lat", "long")) #convert to sf df tibble
class(palestine_cities_sf) #check if class has changed to sf dataframe tibble
## [1] "sf" "tbl_df" "tbl" "data.frame"
st_crs(palestine_cities_sf) #check for geographic projections
## Coordinate Reference System: NA
st_crs(palestine_cities_sf) <- 4326 #assign geographic projection (4326, i.e., longlat WGS84)
st_crs(palestine_cities_sf) #check if the geographic projection was assigned
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
st_crs(palestine_cities_sf)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
palestine_cities_sf #check if "lat" and "long" variables were changed to "geometry" variable.
## Simple feature collection with 244 features and 4 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 31.27 ymin: 34.25 xmax: 32.55 ymax: 35.45
## Geodetic CRS: WGS 84
## # A tibble: 244 × 5
## name country.etc pop capital geometry
## * <chr> <chr> <int> <int> <POINT [°]>
## 1 'Abasan al-Jadidah Palestine 5629 0 (31.31 34.34)
## 2 'Abasan al-Kabirah Palestine 18999 0 (31.32 34.35)
## 3 'Abud Palestine 2456 0 (32.03 35.07)
## 4 'Abwein Palestine 3434 0 (32.03 35.2)
## 5 'Ajjah Palestine 5143 0 (32.37 35.2)
## 6 'Anabta Palestine 7311 0 (32.31 35.12)
## 7 'Anata Palestine 9149 0 (31.81 35.28)
## 8 'Anin Palestine 3717 0 (32.5 35.17)
## 9 'Anzah Palestine 2006 0 (32.37 35.22)
## 10 'Aqbah Jabbar Palestine 6337 0 (31.83 35.45)
## # ℹ 234 more rows
st_write(palestine_cities_sf, "D:/R_Gaza/pal_map/data/palestinian_cities",
driver = "ESRI Shapefile", delete_layer = TRUE) # save the sf data frame as shapefile
## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver
## Deleting layer `palestinian_cities' using driver `ESRI Shapefile'
## Writing layer `palestinian_cities' to data source
## `D:/R_Gaza/pal_map/data/palestinian_cities' using driver `ESRI Shapefile'
## Writing 244 features with 4 fields and geometry type Point.
unzip(zipfile = "D:/R_Gaza/gaza/data/localities_palestine_sf.zip" , exdir ="data/localities_palestine")
localities_palestine_sf <- st_read("data/localities_palestine")
## Multiple layers are present in data source D:\R_Gaza\pal_map\data\localities_palestine, reading layer `%D8%AA%D8%AC%D9%85%D8%B9%D8%A7%D8%AA_%D9%81%D9%84%D8%B3%D8%B7%D9%8A%D9%86___Localities_Palestine'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## automatically selected the first layer in a data source containing more than
## one.
## Reading layer `%D8%AA%D8%AC%D9%85%D8%B9%D8%A7%D8%AA_%D9%81%D9%84%D8%B3%D8%B7%D9%8A%D9%86___Localities_Palestine' from data source `D:\R_Gaza\pal_map\data\localities_palestine' using driver `ESRI Shapefile'
## Simple feature collection with 619 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 3809204 ymin: 3661924 xmax: 3960023 ymax: 3836009
## Projected CRS: WGS 84 / Pseudo-Mercator
localities_palestine_sf
## Simple feature collection with 619 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 3809204 ymin: 3661924 xmax: 3960023 ymax: 3836009
## Projected CRS: WGS 84 / Pseudo-Mercator
## First 10 features:
## OBJECTID REGIONCODE DISTCODE GOVCODE LOCCODE NAMEAR
## 1 1 2 24 2470 24703430 عبسان الجديدة
## 2 2 2 24 2470 24703470 خزاعة
## 3 3 2 24 2465 24653200 مخيم دير البلح
## 4 4 2 24 2465 24653275 وادي السلقا
## 5 5 2 24 2455 24552681 أم النصر (القرية البدوية)
## 6 6 2 24 2465 24653065 مخيم النصيرات
## 7 7 2 24 2470 24703420 خانيونس
## 8 8 2 24 2470 24703425 بني سهيلا
## 9 9 2 24 2470 24703445 عبسان الكبيرة
## 10 10 2 24 2470 24703485 الفخاري
## NAMEEN
## 1 'Abasan al Jadida
## 2 Khuza'a
## 3 Deir al Balah Camp
## 4 Wadi as Salqa
## 5 Um Al-Nnaser (Al Qaraya al Badawiya)
## 6 An Nuseirat Camp
## 7 Khan Yunis
## 8 Bani Suheila
## 9 'Abasan al Kabira
## 10 Al Fukhkhari
## NOTES
## 1 <NA>
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## 7 تشمل المَواصِي (خانيونس)، وقِيزَان النَجَّار، وقَاعْ الخَرَابَة، وقَاعْ القُرِين، وأُمُّ كَمِيل، وأُّمُّ الكِلاب
## 8 <NA>
## 9 <NA>
## 10 <NA>
## Shape__Are Shape__Len geometry
## 1 4919436.7 11010.191 MULTIPOLYGON (((3824072 367...
## 2 9193057.5 14360.481 MULTIPOLYGON (((3826180 367...
## 3 337604.1 3122.714 MULTIPOLYGON (((3823061 368...
## 4 8752155.1 12611.773 MULTIPOLYGON (((3827289 368...
## 5 4306257.9 9547.706 MULTIPOLYGON (((3844490 370...
## 6 1099596.7 10312.221 MULTIPOLYGON (((3827522 369...
## 7 72893548.1 53756.789 MULTIPOLYGON (((3821981 367...
## 8 9181044.3 16897.216 MULTIPOLYGON (((3823716 367...
## 9 17770337.9 20809.409 MULTIPOLYGON (((3823887 367...
## 10 12900017.9 17496.895 MULTIPOLYGON (((3820790 367...
class(localities_palestine_sf)
## [1] "sf" "data.frame"
st_crs(localities_palestine_sf)$proj4string
## [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"
localities_palestine_WGS84 <- st_transform(localities_palestine_sf, crs = 4326)
st_crs(localities_palestine_WGS84)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
st_write(localities_palestine_WGS84, "D:/R_Gaza/pal_map/data/localities_palestine",
driver = "ESRI Shapefile", delete_layer = TRUE)
## Deleting layer `localities_palestine' using driver `ESRI Shapefile'
## Writing layer `localities_palestine' to data source
## `D:/R_Gaza/pal_map/data/localities_palestine' using driver `ESRI Shapefile'
## Writing 619 features with 10 fields and geometry type Multi Polygon.
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: One or several characters
## couldn't be converted correctly from UTF-8 to ISO-8859-1. This warning will
## not be emitted anymore.
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 147785534.120117009 of
## field Shape__Are of feature 57 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 263963660.09277299 of
## field Shape__Are of feature 83 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 135039999.595703006 of
## field Shape__Are of feature 232 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 249401566.34472701 of
## field Shape__Are of feature 323 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 248853481.261718988 of
## field Shape__Are of feature 366 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 300296232.731445014 of
## field Shape__Are of feature 495 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 140128653.509766012 of
## field Shape__Are of feature 497 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 100428703.091796994 of
## field Shape__Are of feature 548 not successfully written. Possibly due to too
## larger number with respect to field width
st_crs(localities_palestine_WGS84)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
st_crs(palestine_cities_sf)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
